Wave chaos in the non-equilibrium dynamics of the Gross-Pitaevskii equation 
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The Gross-Pitaevskii equation (GPE) plays an important role in the description of Bose-Einstein 
condensates (BECs) at the mean-field level. The GPE belongs to the class of non-linear Schrodinger 
equations which are known to feature dynamical instability and collapse for attractive non-linear 
interactions. We show that the GPE with repulsive non-linear interactions typical for BECs features 
chaotic wave dynamics. We find positive Lyapunov exponents for BECs expanding in periodic and 
aperiodic smooth external potentials as well as disorder potentials. Our analysis demonstrates 
that wave chaos characterized by the exponential divergence of nearby initial wavefunctions is to 
be distinguished from the notion of non-integrability of non-linear wave equations. We discuss the 
implications of these observations for the limits of applicability of the GPE, the problem of Anderson 
localization, and the properties of the underlying many-body dynamics. 
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I. INTRODUCTION 

Following the experimental realization of Bose- 
Einstein condensates (BECs) in dilute ultracold gases, 
the Gross-Pitaevskii equation (GPE), has taken cen- 
ter stage to describe the equilibrium as well as non- 
equilibrium dynamics of the condensate at the mean-field 
level 1 . The replacement of the many-body wavefunction 
by the effective single-particle condensate wavefunction 
has proven to be a remarkably successful approximation 
for predicting a large variety of physical observables." 
Among the observables are both ground-state proper- 
ties and elementary excitations in inhomogeneous back- 
ground potentials. 2-9 The GPE belongs to the class of 
non-linear Schrodinger equations (NLSEs) which have a 
broad range of applications ranging from nonlinear optics 
to plasma physics and Bose-Einstein condensation. 10 Ef- 
fects beyond the GPE have been observed in BECs most 
notably in optical lattices with deep wells and small oc- 
cupation numbers per site. In this regime, explicit many- 
body descriptions such as the Bose-Hubbard model are 
more suitable. 11 ' 12 

The non-equilibrium dynamics of BECs, specifically their 
expansion in disordered potentials has recently received 
a lot of attention (see e. g. Ref. 3-9,13-21 and references 
therein). One focus is on the observation of Anderson 
localization of a quantum gas. For weak disorder poten- 
tials (V) <C /i, where (V) is the variance of the potential 
and n is the chemical potential of the BEC, the GPE 
was assumed to be valid during the non-equilibrium ex- 
pansion starting from the BEC released from the trap to 
the dilute localized state for which the linear one-particle 
Schrodinger limit is reached. 15 ' 1 ' The consequences of the 
presence of the non-linearity for the non-equilibrium dy- 
namics in a disordered potential described by the GPE 
deserves a careful analysis. The non-linear Schrodinger 



equation with attractive interactions is known to fea- 
ture dynamical instabilities leading to collapse of the 
wavepacket. 2 Closely related, the GPE with repulsive 
pair interaction in a strictly periodic potential features 
near the Brillouin zone boundary a dynamical (modula- 
tion) instability since the effective negative-mass disper- 
sion translates into an effective attractive pair interaction 
(see e. g. Ref. 22-26 and references therein). Discretized 
models resembling the Fermi-Pasta-Ulam-Tsingou 61 sys- 
tem of non-linearly coupled oscillators have been found 
to feature stochastic dynamics and relaxation with an in- 
creasing entropy. 27-29 

We show in the following that the GPE for realistic pa- 
rameters for the expansion of BECs in the quasi-one- 
dimensional (ID) regime displays true wave chaos as 
measured by a positive Lyapunov exponent in Hilbert 
space. By careful checks of the accuracy of the propa- 
gation including the method of time-reversed propaga- 
tion, this "physical" chaos can be distinguished from the 
numerical chaos previously observed for the NLSE. 27 ' 30 
We furthermore show that chaos goes beyond the non- 
integrability of non- linear wave equations. The physi- 
cal consequences of deterministic chaos in the GPE for 
smooth periodic and aperiodic potentials as well as dis- 
order potentials will be discussed. We argue that wave 
chaos in the GPE is a signature for the breakdown of 
mean-field theory and delimits the border of its applica- 
bility. The latter does not preclude that certain ensem- 
ble expectation values of a BEC can be approximately 
accounted for by a GPE. We conjecture that the chaotic 
fluctuations are a signature of excitations and depletion 
of the condensate. Although our physical interpretations 
focus on BECs, our findings are relevant for other areas 
of application of the NLSE as well. 31 
The paper is organized as follows. In Sec. II we briefly 
describe the model for the expansion of a quasi-lD BEC 
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within the GPE. Numerical methods for the propaga- 
tion of the condensate wavefunctions will be reviewed in 
Sec. III. Non-integrability and wave chaos will be dis- 
cussed in Sec. IV followed by numerical results in Sec. V. 
Conclusions and conjectures will be given in Sec. VI. 



II. GROSS-PITAEVSKII EQUATION FOR AN 
EXPANDING (QUASI)-ID BEC 

The GPE for the inhomogeneous condensate wavefunc- 
tion ip(r, t) — (^(r, i)), the non-vanishing expectation 
value of the field operator ip(r,t) that becomes finite at 
the Bosc-Einstcin phase transition, is given by 



ih 



dt 



fi 2 v 2 

2m 



+ V(r)+g w \iP(r,t)\ 2 ^(r 7 t) 



The effective inter-particle non-linear coupling constant 
in three dimensions, 



ff3D 



47r/i 2 a s 
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is expected to account for the condensate dynamics un- 
der the conditions of low excitations, i.e. weak deple- 
tion of the condensate, and of weak coupling. More- 
over, Eq. 1 assumes short-range interactions at low en- 
ergies such that the particle-particle interaction can be 
described by a contact interaction whose strength is pro- 
portional to the s-wave scattering length a s . The mean- 
field approximation is assumed to be valid in the dilute 
regime n 1 / 3 a s <C 1, where n is the particle density. We 
consider in the following a (quasi)-lD system, for exam- 
ple, a cigar-shaped trap where the radial and longitudi- 
nal frequencies are related as io r 3> oji and we assume the 
chemical potential /i to be small compared to the trans- 
verse quantization energy /i <C hui r - Thus, the BEC is 
described by a ID order parameter. In the transverse di- 
rection the dynamics is confined to the groundstate. The 
ID GPE is then given by 



di(j(x, t) 
1 dt 



2m dx 2 



tp(x,t) + V(x)ip(x,t) 



+2huj r a s N\ip(x, t)\ 2 ip(x, t). 



(3) 



The number of atoms N enters explicitly because we 
normalize the order parameter (in the following termed 
wavefunction) as J dx|i/>(x, t)| 2 = 1. The potential V(x) 
is an external potential to be described in more detail 
below. Measuring energies in units of the longitudinal os- 
cillator energy hui , length in units of the oscillator length 
l = {h/muji) 1 / 2 and time in units t — 1/iJt the GPE 
takes in these (oscillator) units the form 



.di/)(x,t) _ ld 2 *p(x,t) 



dt 



dx 2 



+ U(x)ip(x,t) 



-g\fl>(x,t)\ 2 fl>(x,t), 



(4) 



with U(x) = V(xlo)/h(jJi, g — 2{w T /u){)aN and a = a s /lo- 
In the following we refer to the effective Hamiltonian of 
the linear Schrodinger equation as Hl = ~\~§^2 + U[x) 
and to the non-linear part of the Hamiltonian as -Hnl = 

Most of the numerical simulations presented in the fol- 
lowing are performed for ultracold gases of Rb 87 initially 
stored in a cigar-shaped trap with the following param- 
eters: N = 1.2 x 10 4 and a s = 5.82nm, and following 
Ref. 7 wj = 5.4 x 2n Hz and w r = 70 x 2n Hz. Accord- 
ingly, the unit of time is to = 29.47ms and the unit of 
length is Iq = 4.64/iin. The associated non-linearity is 



g = 2^aN 



390 



(5) 



for later reference. Note that the rather high numerical 
value for g is due to the explicit inclusion of the number 
of particles N and does not contradict the assumption 
of weak coupling. For these parameters /j, tko r is not 
fulfilled, however, it is assumed that the BEC is in the 
quasi-lD regime such that a ID order parameter is suf- 
ficient to describe the dynamics (see e.g. Ref. 6,7,17,18 
where the dynamics is essentially confined to its ground- 
state in the transverse direction). 

The GPE (Eq. 4) can be derived as the Euler-Lagrange 
equation from a Lagrangian functional in the fields 
tp(x, t) and ip(x, t). The corresponding Hamiltonian func- 
tional is given by 



H[i)] = J dx 



\d x i,\ 2 + u\i/jf 
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(6) 



The total energy functional E = H[tp] must be conserved 
during the expansion of the BEC. This conservation law 
provides a stringent test for numerical stability of the 
long-time propagation. 



III. NUMERICAL METHODS 

The propagation of an initial condensate wavefunction 
ip(x,t — 0) according to the GPE (Eq. 4) proceeds by 
discretization of space and time. The space discretiza- 
tion must be performed carefully since simple discretiza- 
tion schemes may unintentionally convert the integrable 
continuous GPE into a discrete non-linear system that 
displays chaos. For example, using the simplest finite 
difference scheme with 5x — Xi + i — x u Eq. 4 with U = 0, 
takes the form 



id t ipi = - 



2^ +1 



2V> l + V I -i)+5l^| 2 V' l - (7) 



Eq. 7 shows chaos 2 ' and has been used to study stochas- 
tic dynamics and thermalization in the mean-field Bosc- 
Hubbard system 29 while its continuous limit is known to 
be integrable. 2 Since we want to study the continuous 
system we have to avoid such discretization artifacts by 
employing a more elaborate discretization. Our spatial 
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discretization is based on the finite element discrete vari- 
able representation (FEDVR) (see Ref. 32-35 and refer- 
ences therein). We split the space into finite elements 
which arc discretized with the help of a discrete variable 
representation (DVR) basis. The basis consists of La- 
grange interpolating polynomials determined via a grid 
consisting of the zeros of Legendre polynomials and the 
endpoints of each element. The elements are connected 
via bridge functions, thus guaranteeing the continuity of 
the wavefunctions. Integrals are approximated via the 
Gauss-Lobatto quadrature. In the FEDVR the local op- 
erators (e. g. the potential operators) are diagonal. The 
kinetic operator within a single element is a full Nb x N\, 
matrix when is the number of basis functions. Within 
the FEDVR the Hamiltonian is a sparse matrix when the 
number of elements and basis functions are chosen appro- 
priately. 

For temporal propagation we tested several different 
propagators: the second-order difference 36 ' 37 (SOD), 
Runge-Kutta 38 , the real space product split operator 37 ' 39 
(equivalent to the symplectic operator 40 SABAi), Crank- 
Nicholson 36 and Lanczos propagators. 3 '' 41 The Runge- 
Kutta and SOD algorithms have proven to be most ac- 
curate in terms of energy conservation. In most of our 
calculations we use the SOD algorithm. It is an explicit 
conditionally stable integration scheme given by the re- 
cursion relation: 

ip(x, t + dt) = ip(x, t-dt)- i2dtHip(x, t), (8) 

where H is given by 

H^-~ + U(x)+g\^x,t)\ 2 . (9) 

The numerically more expensive 4 th order Runge-Kutta 
algorithm is used for cross-checking. Comparison with 
a 5 th order Runge-Kutta algorithm allows an accuracy 
estimate in 2 to be of the order of 10 -15 for each time 
step. 

Both SOD and Runge-Kutta are not symplectic. The 
Hamiltonian structure of the GPE would suggest to use 
a symplectic integrator. Such integrators preserve, by 
design, the volume in phase space which is characteristic 
for the dynamics of Hamiltonian systems. However, it 
has been demonstrated for the NLSE and the Korteweg- 
deVries equation and conjectured for infinite dimensional 
Hamiltonian systems in general that the use of symplec- 
tic time integrators is less important than high accuracy 
of spatial derivatives. 42 ' 43 For the results presented in the 
following we use the time-increment dt — 4 x 1CP 6 and 
5 basis functions for finite-element widths of typically 
dx = 0.08 in a box of size x G [-1000, 1000]. (For the 
long time propagations presented in Sec. V C the box size 
is x 6 [—1600, 1600].) We use hard wall boundary condi- 
tions for our numerical box. The numerical box is cho- 
sen sufficiently large such that the spreading wavepacket 
ijj{x,t) does not effectively reach the walls of the box 
within the time of propagation. 



Apart from energy conservation, an additional sensitive 
test of the numerical stability of the propagation algo- 
rithm is the time reversal of propagation. Expressing 
the propagation of the initial condensate wavefunction 

ip(x, 0) as 

1>(x,t) = U t [il>(x,Q)], (10) 

the time reversed propagation 

xP(x,0) = U^(x,t)] (11) 

should recover the initial state. For the linear 
Schrodinger equation (LSE) , reaching the initial state via 
Eq. 11 is easier than for the NLSE due to the intrinsic sta- 
bility of the linear dynamical evolution. For the NLSE, 
Eq. 11 provides a measure for the accumulation of nu- 
merical noise, i. e. a measure for "numerical chaos" . 30 
In order to distinguish true deterministic ("physical") 
chaos from numerical noise, we will probe the accuracy 
of Eq. 11 by quantifying how close the backward propa- 
gated wavefunction will be to the initial state if)(x, 0). 
Another stringent test results fom the scaling property 
of the GPE. We subject the space and time coordinates 
to the scaling (7 > 0) 

t-ti = t/7 

x ->• x = z/Vt- (I 2 ) 

The GPE becomes, after multiplication by 7, 

d 1 d 2 _ -__ 

i-^ip(x,t) = --■^ip(x,T) + U(x)i/)(x,t) 

+ ig \^(x,F)\ 2 ^(x,t), (13) 

with U(x) = jUl^y^yx). The normalization of the con- 
densate wavefunction after rescaling 

#g,i)=7 V V0M~) (14) 
leads to the rescaled GPE 

d - _ l d 2 - _ 

i-^tp(x,i) = ---^ip(x,i) + U(x)?p(x,t) 

+ gfr)$(x, t)\H{x,t). (15) 

with (7(7) = y/jg. Converged numerical solutions that 
are not subject to numerical noise will satisfy the scaling 
behavior (Eq. 15) for every positive value of 7. Scal- 
ing is equivalent to a continuous variation of the spatio- 
temporal grid. One further consequence of Eq. 15 is that 
scaling corresponds to free tuning of the non-linearity 
in the absence of U. Therefore, thresholds for critical 
strength of non-linearity, if they exist at all, can occur 
only in the presence of external potentials. 

IV. INTEGRABILITY, NON-INTEGRABILITY 
AND WAVE CHAOS IN THE GPE 

For non-linear wave equations, the GPE being an ex- 
ample of which, integrability is closely associated with 
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the existence of solutions via the inverse scattering trans- 
form (1ST). Briefly, the initial value problem (IVP) of 
the non-linear wave equation possesses an integrable so- 
lution if the solution can be determined by 1ST. In such 
a case, the non-linear wave equation is associated with 
a system of auxiliary linear ordinary differential equa- 
tions (LODEs) for which ip(x,t) acts as a "potential" 
and whose direct and inverse scattering solution solves 
the original IVP (see e. g. Ref. 44-46). The time t acts 
here as a "deformation parameter" . While this criterion 
for integrability is explicit and constructive, it does not, 
however, provide a priori a sufficient criterion for non- 
integrability. 

For the GPE, a few special cases are known: For van- 
ishing external potential, U(x) = 0, and infinite do- 
main (—oo < x < oo), the GPE is integrable and 
the corresponding LODEs are the Zakharov-Shabat (ZS) 



system. 



While for attractive inter-particle interac- 



tions (g < 0), the ZS system features both a continuum 
and bound states corresponding to soliton solutions of the 
GPE, for repulsive interactions (g > 0), the case of inter- 
est in the following, the ZS spectrum is purely continu- 
ous corresponding to a dispersing (decaying) condensate 
wavefunction ij){x,t). For hard- wall Dirichlet boundary 
conditions, the GPE with U{x) — was conjectured to 
be non-integrable. 28 

For non- vanishing external potentials (U(x) ^ 0), no 
general results on the existence of integrable solutions 
of the GPE are known. For a few special cases, inte- 
grability has been demonstrated. These cases include 
the harmonic potential, U(x) cx x 2 , with time-dependent 
non-linearities, 47-49 and linear potentials, U(x) cx x, with 
time- independent non-linearity. 50 ' 51 
One test for integrability is the Weiss- Tabor-Carnevale 
(WTC) test for non-autonomous NLSEs with time 
and space dependent dispersion, non-linearity and 
dissipation. 45 The WTC test is based on the Painleve 
conjecture (see e.g. Ref. 45,46) and provides a necessary 
condition for integrability (see Ref. 52 for specific poten- 
tials and nonlinearities) . Another integrability condition 
is based on the Lax pair method yielding, in general, dif- 
ferent criteria. 03 However, both integrability conditions 
agree in that for constant dispersion and nonlincarity and 
vanishing dissipation the potential U(x) may be at most 
linear in x. For disordered potentials in the NLSE it 
was shown that quasi-periodic motion may persist for 
weak non-linearities that can be considered as a small 
perturbation. 5 

In the following we consider potentials U (x) of physical 
interest which do not satisfy the known criteria of integra- 
bility. They include disordered, periodic and aperiodic 
potentials. We will characterize the resulting properties 
of the GPE by a stronger criterion than non-integrability, 
i. e. the appearance of wave chaos. As a measure for 
wave chaos we employ the existence of a positive Lya- 
punov exponent which indicates exponential sensitivity 
to initial conditions. In analogy to classical (particle) 
chaos, where the distance of initially close trajectories 



is followed in phase space, we follow the distance of ini- 
tially close wavefunctions in Hilbert space. The metric 
in Hilbert space is the L 2 norm. Accordingly, we use a 
distance function 

d«tyi,lM) = ^i-Wi-^2> (16) 



dx \ipi(x,t) - ip 2 (x,t)^ 



The Lyapunov exponent follows as the limit 
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A = — lim 



lim 



In 



(17) 



In analogy to chaos for point particles, A is only well- 
defined in terms of a coupled double limit. The distance 
function (Eq. 16) for normalized functions and norm- 
conserving evolution reduces to 



d (2) (V>i,^ 2 ) = l-Re(Vi|^ : 



2 ■ 



(18) 



It is bounded by < df® < 2. The upper bound corre- 
sponds to ipi = —tp2- The value 1 is reached for complete 
orthogonality and thus corresponds to the "maximal" 
separation in Hilbert space. Therefore, we expect that an 
exponential separation of initially "nearby" , i. e. almost 
identical wavefunctions will eventually saturate at values 
e?( 2 ) 1. Eq. 17 can be viewed as a measure for the expo- 
nential separation of trajectories in Hilbert space on the 
hypersphere S of unit radius defined by wavefunctions 
of unit norm. In practice, A is numerically determined 
by the slope of that segment of growth on a semi-log 
plot that is approximately linear. The choice of initial 
displacements d^ 2 '(ipi, ip2',t) admitted in Eq. 16 is con- 
strained by the Hamiltonian structure of the GPE with 
its corresponding energy shells, i. e. H[ipi] = H [^2] = E. 
With this constraint, the evolution of ipi, ip2 proceeds on 
a hypersphere in Hilbert space of fixed E. 
It is instructive to first explore the short-time behavior 
of d( 2 \ipi,ip2',t). Point of reference is the LSE, for which 
d (2 ^ is strictly conserved, d^ 2 \ij)i^2]t) — d^(ipi,ipi]0). 
Variation of d^ is therefore a direct measure for the in- 
fluence of the non-linearity. For the non-linear GPE, we 
obtain for initial real wavefunctions ipi(x, 0) and ip2{x, 0) 
up to second order in t 



Re(Mt)\Mt)) = c + c 2 t 2 + 0(t 4 ) 



(19) 



with 



and 



c = (Vi (0)|V 2 (0)> = l-d (2) (^i,feO) (20) 



c 2 = 2 / dir V'l( a; ) Q )[«l( af :0) - Tl2(x,0),Hj J ]lp 2 (X!,0) 



^ / dx ipi(x,0)( ni (x,0) - n 2 (x,0)yip2(x,0) 



(21) 
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FIG. 1: (Color online) Distance function d^ 2 \^i,tp2',t) 
(Eq. 16) for a freely expanding BEC (U(x) = 0) with non- 
linearity go released from a harmonic trap (blue solid line). 
The initial normalized states tpi, ip2 are slightly perturbed 
ground states of the interacting system in the harmonic trap 
with a small linear perturbation (same states as will be used 
and explained in more detail in Eq. 23 in Sec. V). The initial 
distance function is d (2) (0) w 2.8 x 10~ r . The black dashed 
line corresponds to w- 2 ' (ipi^^t) according to Eq. 18 and 
Eq. 19. 

with m(x,0) = \ipi(x, 0)| 2 the condensate density. The 
square bracket in the first term of Eq. 21 denotes the 
usual quantum mechanical commutator. For strongly 
interacting condensates where the interaction energy is 
large compared to the single-particle energies, the g 2 
term in Eq. 21 dominates. The initial quadratic increase 
of the distance function S- 2 ' with time is universal and 
independent of integrability, non-integrability, or wave 
chaos. The short-time behavior of the distance function 
for a freely expanding condensate (U(x) = 0) released 
from a harmonic trap displays the growth with t 2 fol- 
lowed by saturation (Fig. 1). We conjecture that this 
saturation is a signature of the integrability of the GPE 
for U(x) =0. 

A slightly different problem, expansion of a conden- 
sate in zero potential (U(x) = 0) inside a box de- 
fined by Dirichlet boundary conditions, was investigated 
by Villain and Lewenstein. 28 Motivated by the struc- 
tural similarity to the famous Fermi- Pasta- Ulam-Tsingou 
model, they observed stochastic behavior as extracted 
from a normalized statistical entropy. For sufficiently 
strong non-linearity a trend toward equipartition of en- 
ergy among the modes of the LSE was found. The au- 
thors identified the free GPE with Dirichlet boundary 
conditions as non-integrable. We have investigated the 
behavior of the distance function for this system. 
We follow Ref. 28 and chose as an initial state a Gaus- 
sian wavepacket iI)g(x) with width a = 1 in a box of 
length L = 20. This initial state is slightly distorted by 
a small linear admixture of strength a: 



with normalization constants N\$- We test the sensi- 
tivity of d^(t) to a [i.e. to d^ 2 '(0)] by choosing two 
different values for a: a = 5 x 10 which gives a sim- 
ilar initial d (2) (0) as in Fig. 1 and a = 5 x 10 5 for 
which g?( 2 )(0) is two orders of magnitude smaller. The 
numerical results confirm the expectation that d^ 2 \t) is 
linearly proportional to the initial distance d^ 2 \0) (con- 
stant shift in a logarithmic plot, see Fig. 2). For the 
strength of the non-linearity we consider both a moder- 
ately strong value g/L ss 79ei where t\ — ir 2 /2L 2 is the 
ground-state energy of the LSE in the box comparable to 
the value (g/L 61ei) in Ref. 28, and a much stronger 
non-linearity (g/L w 1582ei) corresponding to go (Eq. 5). 
The energy unit in Ref. 28 is related to our energy unit 
as e 1 = 0.0123/k^. The time unit T x = 4L 2 /n « 500t is 
long compared to our "oscillation" time scale to charac- 
terizing the initial state. 

For both non-linearities the distance function d^ (t) in- 
creases quadratically with a slope approximately propor- 
tional to the non-linearity (see Fig. 2) in qualitative ac- 
cord to the short-time behavior (Eq. 18 and Eq. 19), how- 
ever with much reduced slope. The growth seen in Fig. 2 
proceeds slowly. The quadratic increase is fundamentally 
different from the exponential increase observed in the 
following sections. This allows identifying within non- 
integrable systems those that do not exhibit wave chaos, 
as in Ref. 28, and those that do, as the systems discussed 
in the following. 

The overall quadratic rise is modulated by fluctuations 
which result from partial revivals of the underlying LSE 
(see inset of Fig. 2). One prominent oscillation period 
is Ti/8 corresponding to the inverse energy spacing be- 
tween the two lowest LSE states of even parity. 28 Another 
prominent period T\ /4 results from the classical periodic 
orbit in the box. The non-linearity leads to a prolifera- 
tion of frequencies by side-band coupling. 
The important conclusion is that the non-integrable dy- 
namics of the free GPE in a box with hard walls is, 
despite its complex appearance and the trend toward 
equipartition among the linear modes, non-chaotic as 
measured by the exponential separation of nearby tra- 
jectories in Hilbert space. Chaos requires the presence of 
non- vanishing external potentials U(x). 

V. NUMERICAL EXAMPLES FOR WAVE 
CHAOS 

We present in this section examples for expansion of 
BECs in potentials U(x) of experimental relevance. For 
the generation of nearby initial states we use the following 
scenario: The BEC is created in a harmonic trap, i. e. the 
condensate wavefunction is given by the groundstate of 
the GPE, ipg( x ), which is close to the Thomas- Fermi pro- 
file for the parameters chosen (see Sec. II). Two nearby 
normalized states are then created by weak perturbations 
linear in x, 



iM&.O) =JVi, 2 (l±oa;)^G(aO (22) 



lMa:,0) =N lj2 (l±ax)i) & (x), (23) 
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with periodicity I and strength of the weak potential Uq. 
We require 



t/Ti 



FIG. 2: (Color online) (a) Log-log plot of the distance func- 
tion 

d( 2 >(Vi,V2;t) (Eq. 16) as a function of t for two initially 
nearby Gaussian wavepackets in the box of length L obtained 
from propagation of the GPE with moderate non-linearity 
(g/L ~ 79ei) and strong non-linearity (g/L m 1582ei). For 
g/L « 1582ei two different distortion parameters a (see 
Eq. 22) were used, a — 5 x 10 -4 and a = 5 x 10 -5 (marked 
by arrows in the figure). Time is measured in units of 
T% = 27r/ei. For reference, the exponential growth in case 
of wave chaos with a Lyapunov exponent as in Fig. 4 is 
shown as dotted line, (b) The curve with g/L w 1582ei and 
a = 5 x 10~ 4 in linear scale. 



with normalization constants N\^- Simultaneously with 
the release at t = the external potential U(x) is 
switched on. For the perturbation parameter a control- 
ling the distance between ipi and ip2 we typically choose 
a = 10~ 4 . We have verified that the extracted values for 
the Lyapunov exponent A are independent of the value 
of a. 



A. Weak periodic potential 

We first consider the expansion of the BEC in weak 
periodic potentials 



U /e « 1, 



(25) 



where e = E/N is the total energy per particle of the in- 
teracting system. We will use in the following Uq = 0.2e. 
The chemical potential is to a good degree of accuracy 
fi(t = 0) = 2e at t = for our choice of initial conditions. 
Unlike e, /i(£) is not conserved during the expansion and 
therefore, we use e as characteristic energy scale. The pe- 
riodic potential gives rise to a band structure of the linear 
problem with bandwidth of the first band, W = ir 2 /2l 2 , 
(in literature also called recoil energy 26 ) and a band gap 
at the first Brillouin zone boundary, A ss Uq- The wide 
band gap limit corresponds to W/U -C 1, the nearly free 
particle limit to W/Uo 3> 1. The spatial periodicity I is 
assumed to exceed the coherence (healing) length £ of the 
condensate 



1>Z 



with 



8c 



(26) 



(27) 



U(x) = Uq cos (2nx/l) 



(24) 



In the opposite limit £ <C I, the condensate can no longer 
resolve the local variation of the potential. 
The initial wavefunctions are chosen as above (Eq. 23) 
having identical energy in the potential (Eq. 24) and a 
Thomas-Fermi length of Ltf ~ 8.4 which is larger than 
the potential periodicity I = 5.8£ with £ = 0.0945. The 
non-linearity is go f=s 390. For t > the condensate 
wavefunctions expand and acquire an increasingly com- 
plex fluctuation pattern. Appreciable deviations start to 
emerge near the center x = where fluctuations have 
the largest amplitude [see Fig. 3 (a)]. The deviations 
increase with increasing time [see Fig. 3 (b)], such that 
locally the amplitudes of and tp2 become uncorrelated 
[see the inset of Fig. 3 (b)]. Such an extreme sensitivity 
to initial conditions is the hallmark of wave chaos. It 
is worth noting that on the length scale of / both wave- 
functions tend to mimic the oscillations of the potential 
while on larger length scales fluctuations of ip\ and ?/>2 
lose their correlation. 

We quantify the loss of correlation and the seemingly ran- 
dom fluctuations by the Lyapunov exponent A (Eq. 17). 
Positive A signifies wave chaos. The time evolution of 
the distance function (ipi , V>2! t) resembles the time 
evolution of particle chaos. After an initial generic non- 
exponential (quadratic) growth (see Fig. 4), dS 2 ' grows 
exponentially until the regime of fluctuations around the 
saturation value is reached. The slope of the exponen- 
tial growth allows to numerically extract the Lyapunov 
exponent (A ps 0.7 for the case of Fig. 4). 

The deterministic nature of the the exponential sepa- 
ration can be verified by the time reversal test (Eq. 11). 
Upon reversal of the direction of time the two conden- 
sate "trajectories" approach each other again to distances 
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FIG. 3: (Color online) Onset of divergence between ipi and 
ip2 as a function of time. ^i] 2 and \ip2\ 2 are given for (a) 
t = 15 and (b) t = 20. Inset in (a): comparison of the 
spatial fluctuations of the wavefunction and the period of the 
potential U(x). The potential parameters are Uq = 0.2e and 
I — 5.8£. The non-linearity is go- Inset in (b): increasingly 
uncorrelated fluctuations of l^il 2 an d |^2| 2 - 



close to their initial value (Fig. 5). The distance function 
between ipi(x, 0) and U_ t [ipi(x, t)] is of the order of 10~ 9 . 
The near-perfect reversibility precludes numerical noise 
or "numerical chaos" as a source for divergence. 
Furthermore, the scaling property of the GPE (Eq. 15) 
is verified to a high degree of accuracy (Fig. 5) . We scale 
the initial conditions according to p,(t = 0) = y(J,(t = 0), 
such that U/e is constant. Likewise the period of the po- 
tential is scaled such that f/£ = y/jl/y/j£, is kept fixed. 
Despite a large absolute variation of S- 2 ^ with 7 [Fig. 5 
(a)], upon rescaling, d/ 2 ^ coalesces to a remarkable degree 
of accuracy [Fig. 5 (b)]. The 7 scaling is also equiva- 
lent to an effective change of the numerical grid in (2, i), 
i. e. an increase of the numerical accuracy, since we have 
used the same absolute grid size for all calculations. The 
invariance of the results clearly demonstrates numerical 




FIG. 4: (Color online) Exponential growth of d' 2 ' for two con- 
densate wavefunctions propagated by the GPE in a weak peri- 
odic potential (same parameters as in Fig. 3). Time-reversed 
propagation for t > to recovers the initial state precluding 
numerical noise as origin of the exponential growth. Two dif- 
ferent distortion parameters have been used: a = 10 -4 and 
a = 1CT 5 . 



stability. Since the Lyapunov exponent scales inversely 
with time A = 7A, A decreases with decreasing 7 in accor- 
dance with Fig. 5 (a). Note, however, that this power-law 
scaling does not provide information on the existence (or 
absence) of critical thresholds for wave chaos (see below) . 
In the single-band mean-field Bose Hubbard model, 
a threshold for stochastic dynamics has recently been 
observed. 29 We explore now the behavior of A in the 
continuum analogue, the GPE with a periodic potential. 
The Lyapunov exponent A will, in general, depend on the 
nonlinearity g, the period I and the amplitude Uq of the 
potential. The corresponding three energy scales are the 
total energy per particle e (controlled by the non-linearity 
g) the bandwidth of the first band W, and the ampli- 
tude Uq. The scaling property allows to interrelate the 
dependences of A on these three parameters. For exam- 
ple one can investigate the dependence of A on the ratio 
e/W at constant Uq by keeping e constant and decreas- 
ing W (increasing I) giving the Lyapunov exponent A( c ), 
or keeping W constant and increasing g, hence, e giv- 
ing X(wy If the ratio e/W remains invariant for the two 
cases, the Lyapunov exponents are interrelated by the 
scaling property (provided one of the initial wavefunc- 
tions is rescaled) as A( c ) = j\(w) with 7 < 1. Therefore, 
it is sufficient to consider the ratio e/W and not W and 
e separately. We arrive at a two-dimensional parameter 
plane (Uo/e, e/W). A scan of A is performed from the 
first to the second band of the linear system. 
We observe a clear threshold of A in form of a rapid in- 
crease along both the Uo/e and the e/W axis (see Fig. 6). 
Below the threshold A is vanishing. A power-law fit to 
the threshold line gives e/W cx (U /e)~ - 65 (see dashed 
line in Fig. 6). As the control parameter e/W can be 
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FIG. 5: (Color online) Distance d' 2 ' between Tpi and Tp2 as 
a function of time for different 7. The parameters of the 
periodic potential are U = 0.2e and I = 5.8£. (b) Logarithmic 
plot of (a) squared after rescaling time according to t — yt. 
The curves coincide. The gray solid lines underline the region 
of exponential growth and saturation and are guides for the 
eye. 




FIG. 6: (Color online) The Lyapunov exponent A as a function 
of e/W and Uo/e. e is the energy per particle, W = n 2 /2l 2 
is the bandwidth of the first band of the linear system, and 
Uo is the amplitude of the potential (Eq. 24). The period 
/ is kept fixed / = 0.54811. The same initial wavefunction 
has been used for all data points. A has been determined 
numerically via a linear fit to the logarithmic increase of d (2) 
(Eq.16). The deviation (summed absolute difference between 
the linear fit and the numerical curve) is on average well below 
0.4 and maximal « 0.6. To improve the graphical appearance 
every second point in the Figure is an interpolation between 
two numerical points. The white dashed curve is a fit to the 
threshold according to e/W = 0.165(l7o/e)~ - 65 . The white 
dots mark the point where the energy e equals to the energy 
at the point of inflection of the first band. 



Uo/e. More importantly, the threshold does not seem to 
be sensitive to specific features of the linear problem, as 
presented in the next section. 



B. Weak aperiodic potential 



related to the parameter k, the ratio between the non- 
linearity parameter and the hopping parameter of the 
discrete system, it is suggestive to assume that the tran- 
sition in the discrete system 29 and in the present continu- 
ous system are of common origin. However, the discrete 
system does not account for the transition along Uo/e. 
In other words, the mapping to the single-band Bosc- 
Hubbard is not straightforward and a comparison to a 
multi-band model seems to be necessary. 
Exponential divergence of wavefunctions above a thresh- 
old resembles the dynamical (or modulation) instability 
previously observed. 22 ' 24,26 The modulation instability 
can be related to energetically surpassing the point of in- 
flection along the first band of the linear problem: when 
the second derivative of the dispersion relation turns 
from positive to negative. 24 The position of the inflec- 
tion point, however, does not agree with the position of 
the threshold (in Fig. 6 marked by white dots). The 
reason may be that the simplified picture of modulation 
instability also does not account for the threshold along 



In order to explore whether wave chaos in the GPE is 
a specific feature of the interplay between the periodic 
linear problem and the non-linearity, we investigate the 
dynamics in a strictly aperiodic potential 

U{x) = Uo [ci cos (2-Kx/h) + c 2 cos {2irx/h]\ , (28) 

where li/l\ = 4(1 + is the "most" irrational golden 
mean number and c\ — ci — 1/2. A well-defined band 
gap or well-defined negative mass dispersion relation near 
the Brillouin zone boundary which invoked the explana- 
tion of instability and stochastic motion 22 ' 24 ' 26 are ab- 
sent for Eq. 28. For comparable values of Uo/e as in 
the periodic case we find positive Lyapunov exponents 
(see Fig. 7). Interestingly, the threshold as well as the 
overall behavior agree well with the periodic case, sug- 
gesting that chaos does not depend on the specifics of 
the linear system but rather on the overall length and 
energy scales. In agreement with Fig. 6 the threshold 
shifts toward smaller values of Uo/e for higher ratios of 
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FIG. 7: (Color online) Lyapunov exponent A as a function of 
Uo/e for the periodic and aperiodic potential. The period of 
the periodic potential is (a) I = 0.54811 and (b) I = 0.94503. 
h of the aperiodic potential (Eq. 28) is chosen respectively 
h = I- The potentials are given in the inset. The energy 
e ~ 14 is constant such that e/W is larger in (b) e/W ~ 2.5 
than in (a) e/W « 0.85. 



e/W [Fig. 7 (a) and (b)]. Note that a characteristic en- 
ergy W can be denned for any potential with a charac- 
teristic length l c as W oc l/l 2 . Obviously, the wave chaos 
we have identified in the GPE is a general feature not 
coupled to specific properties of the bandstructure of the 
periodic linear problem. 



C. Weak disorder 

A further generalization is the GPE with random dis- 
order replacing the smooth (a)periodic potential. This 
case is of particular interest in the context of Anderson 
localization recently studied for expanding BECs. 7 ' 8 The 
latter experimental observations refer to the (quasi) lin- 
ear regime when the non-linearity was sufficiently weak 
such that the localization dynamics is governed by the 




FIG. 8: (Color online) Disorder potential for different values 
of the correlation length a. The amplitude for both poten- 
tials is Uo = 0.2e and Nd = 2 x 10 4 . The energy e and healing 
length £ correspond to the groundstate of the GPE in a har- 
monic oscillator at non-linearity go (Eq. 5). 



Schrddinger equation. Localization in the presence of 
interactions has remained an open question. For the dis- 
crete NLSE, subdiffusive expansion rather than localiza- 
tion was observed numerically 20 ' 55,56 We focus on the 
interplay between the randomness induced by the poten- 
tial and by the wave chaos in the GPE. 
We follow closely the scenario employed in the investi- 
gations of Anderson localization: 7 ' 17 the BEC initially 
trapped in a harmonic oscillator expands in a disorder 
potential. We construct a disorder potential with Gaus- 
sian correlation and zero mean. At equidistant grid 
points Xi we place Gaussians with width a and random 
weight Ai 



N d 

1=1 



(29) 



After subtracting the mean, AU(x) — U(x) — (U(x)) and 
normalizing by the standard deviation, (AU(x) 2 ) 1 ^ 2 , the 
disorder potential becomes (Fig. 8) 



U d (x) 



Uo 



(AU(x) 2 ) 1 / 2 



AU{x). 



(30) 



For Uq we choose the same values as in the cases of 
periodic and aperiodic potentials above. The disorder 
potential Eq. 30 is by construction Gaussian correlated, 

C(x) = (U d (x + x)U d (xo)) = U 2 ei^ (31) 
with a Gaussian shaped Fourier transform (Fig. 9) 

C(k) = V2<rU 2 e- k2 ' y2 . (32) 
Unlike for speckle potentials, the averages over odd pow- 
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\C(k)\* 
smoothed 
analytical formula 
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FIG. 9: (Color online) The Fourier transform of the corre- 
lation function C(k) for U = 0.2e and a = 0.7£ (Fig. 8). 
C(k) agrees very well with Eq. 32 when smoothed over the 
fluctuations. 



crs (U^ n+1 ) vanish (numerically only approximately). 
The correlation length of the disorder potential a pro- 
vides a length scale competing with the healing length £ 
of the free condensate. 

For a < £ the wavefunction is exponentially localized 1 ' 
which is characteristic for Anderson localization [see 
Fig. 10 (b)]. For a > £ (we have chosen a = 2£), the 
condensate wavefunction develops algebraically decaying 
tails 17 that travel seemingly undisturbed through the dis- 
order potential similar to propagation in free space [see 
Fig. 10 (a)]. 

We probe for wave chaos for these two cases. It should 
be noted that in the presence of disorder, the determi- 
nation of the exponential separation of nearby initial 
wavefunctions on the energy hypersphere faces the ad- 
ditional difficulty that it is not straightforward to cre- 
ate two nearby wavefunctions tpi(x, 0) and ip2(x, 0) with 
equal energy. In the present case the energies of the wave- 
functions are kept close with a relative deviation of the 
order of 10~ 4 . 

In both cases, the algebraic and the exponential regime, 
we observe an exponentially growing separation and 
eventual saturation near the maximal value of separation 
dP^{ipx,ip2) ~ 1, confirming the intuitive notion that a 
random potential induces wave chaos. The Lyapunov ex- 
ponent is A w 0.85 in the algebraic regime [Fig. 10 (a)] 
and A w 0.86 in the exponential regime [Fig. 10 (b)]. 
It is instructive to inquire into the relation between wave 
chaos and localization. Wave chaos is a signature of the 
non-linearity in the wave equation, while localization is 
the hallmark of disorder in the linear wave equation. In 
discrete non-linear models a variety of relations have been 
observed, one of which is subdiffusive growth rather than 
localization (see e. g. Ref. 56 and references therein). As 
measure for the localization dynamics we use the time 
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FIG. 10: (Color online) Chaotic fluctuations of ip\ and ip2 
at t = 24 for (a) a = 2£ in the algebraic regime and (b) for 
a — 0.7£ in the "localized" regime (in the potentials of Fig. 8). 
The local fluctuations of and ip2 strongly deviate. 



evolution of the variance of the condensate wavefunc- 
tion, Ax — <J (x 2 ) — (x) 2 . Note that our measure for 
the width is the square root of the second moment m-i 
in Ref. 56. We observe (Fig. 11) for the expansion in a 
disorder potential a slowing down of the spread which lo- 
cally scales like Ax cx t a . In qualitative agreement with 
the prediction in Ref. 56 the exponent is not constant 
but changes over time. We observe for the system in 
Fig. 11 with non-linearity go an exponent of a ~ 0.75 up 
to t = 60 (compare our data also with Ref. 17). To stay 
numerically on the safe side during long-time evolutions 
we performed also calculations with go/10. The expo- 
nent is predicted to be independent of non-linearity. 56 
Note that for a direct comparison the time in units used 
in this manuscript is to be scaled by a factor 62 1250 to 
obtain the units of Ref. 56. We observe that the exponent 
a changes from 0.78 between t = 0— 100 to 0.52 between 
t = 100 - 300 and becomes 0.34 between t = 300 - 400 in 
qualitative agreement to Ref. 56. It should be noted, 
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FIG. 11: (Color online) The variance Aa; of the wavefunction 
tpi(x,t) for propagations in free space, in a periodic potential 
(Uo = 0.2e, I — 20£) and in a disorder potential (Co = 0.2e 
and a = 0.7£ as in Fig. 8). The black dashed lines correspond 
to ip2{x, i) for the periodic and disorder potentials where chaos 
is present. The non-linearity is go « 390. The inset gives 
a long time evolution for a weaker non-linearity go/ 10 in a 
disorder potential with Uo = 0.2e and a = 0.32£. In this case 
the initial wavefunction corresponds to the groundstate of the 
GPE for 30/IO. 



however, that the extraction of this exponent for the 
chaotic GPE is numerically much more challenging and 
less accurate than for discrete models. For example, the 
rigorous convergence test for time-reversed propagation 
(Eq. 11) begins to fail for the longest propagation times 
displayed in Fig. 11. 

Without spatial confinement of the system the density 
may become eventually sufficiently low such that the non- 
linearity (gl'iAl 2 ) can be neglected compared to the disor- 
der potential and the system reacts by localization in the 
linear regime. The experimentally observed Anderson lo- 
calization of expanding BECs' is likely associated with 
the approach of the linear regime. However, the presence 
of subdiffusive expansion due to residual non-linearities 
cannot be ruled out. 

A further interesting observation is that the variance 
(Ax) for the two nearby initial states ipi and ip2 agrees 
well (Fig. 11). Despite the intrinsic randomness of the 
propagated condensate wavefunctions on length scales of 
the potential, do averages such as Ax agree very well 
with each other. 



VI. DISCUSSION AND CONCLUSIONS 

We have shown that the expansion of the conden- 
sate wavefunction described by the GPE features wave 
chaos as determined by a positive Lyapunov exponent 
in Hilbert space. We find wave chaos to be generic, 
i. e. present for a variety of one-body potentials U(x) in- 
cluding weak periodic, aperiodic, and random potentials. 
The wavefunctions develop on the length scale of the 
healing length and of the local variation of U(x) fluctua- 
tions with amplitudes and phases exponentially sensitive 
to the initial state such that nearby states become ap- 
proximately orthogonal to each other. We note a remark- 
able exception: for harmonic potentials, U(x) = \x 2 , 
wave chaos is absent as a consequence of the harmonic 
potential theorem for many-body states. 
The physical implications of these findings are of con- 
siderable interest and raise many conceptual questions. 
Wave chaos is a signature of the development of con- 
densate density fluctuations. For a weakly interacting 
Bose gas condensate density fluctuations coincide with 
total density fluctuations of the many-body system and 
represent long-wave phonon excitations. For stronger in- 
teractions, the delayed onset of exponential divergence 
may delimit the characteristic time over which the GPE 
for the condensate is capable of describing the expanding 
BEC. Beyond this time, the expanding condensate is de- 
pleted by multiple excitations which the GPE attempts 
to mimic within a mean-field one-particle wavefunction 
by random fluctuations. Simply put, in this regime wave 
chaos may mark the breakdown of the mean-field ap- 
proximation. Such reasoning would be based on the 
lack of exponential sensitivity to initial conditions in the 
real many-body system. If we consider two many-body 
wavefunctions for N particles, ipi(x±, . . . ,xw,t = 0) and 
4>2(xi, ■ ■ ■ , xn, t — 0), expanded in the basis of energy 
eigenfunctions tp n (xi, . . . ,xn) of the underlying Hamil- 
tonian: ipifei, . . . , x N , t = 0) = J2 n c n<l>n(x\, . . . , x N ) 
and ip 2 (xx,...,x N ,t = 0) = J2 n djpnixi, . . . , x N ) the 
distance function d^^i, i/^) remains constant irrespec- 
tive of the choice of the underlying Hamiltonian. The 
overlap of wavefunctions does not change in time. This 
fundamental property of linear time evolution should be 
reproduced by any approximative theory. However, the 
difficulty with this argument lies in the fact that the 
distance function S 2 ^ in the many-body Hilbert space 
% = YliLiiL 2 ) 1 is unrelated to that in the L 2 space of 
the reduced mean-field wavefunction entering the GPE. 
Moreover, ensemble expectation values, such as the mo- 
ments of the distribution (e. g. Ax discussed above) 
which represent averages over the condensate wavefunc- 
tion over fine-scale density fluctuations, are found to be 
insensitive to small variations of the initial state. The 
GPE may remain predictive for such expectation values 
on longer times scales. The latter would explain the large 
number of successful applications of the GPE to conden- 
sate expansions (see e.g. Ref. 3,7-9,13). Such a link 
between chaotic wavefunctions and observables would 
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be analogous to classical chaos for particles and phase 
space distribution functions: While the long-time evo- 
lution of individual trajectories become unpredictable, 
the ensemble or time average over stochastic regions in 
phase space remain well defined and yield stable ensem- 
ble expectations values. Going beyond large-scale aver- 
aging, the characterization of the fine-scale random fluc- 
tuations developing under propagation with the GPE re- 
mains a widely open problem. Propagation of many- 
body wavefunctions beyond the GPE may shed light on 
the question to which extent the GPE fluctuations, in the 
mean-field level, may faithfully represent at least some of 
the features of multi-particle excitation. The latter may 
also provide insight into the applicability of a wavefunc- 
tion "thermalization hypothesis" originally developed for 
wavefunctions of the LSE for classically chaotic systems 58 
to wave chaos in the GPE. Furthermore, it would be of 
interest to experimentally search for and characterize lo- 
cal fluctuations and structures in expanding condensates, 
e. g. by elastic light scattering. 
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